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Final Report on (Irant NSG5154 

A Study of Autonomous Satell ite N avi sat ion Methods 
Usins the Global Posltionins Satellite System 

Introductio n 

The. development of the NAVSTAR f.lobal Positioning System (GPS) 
will allow satellites to perform orbit determination calculations in 
real time with on-board computers. Special orbit determination algor- 
ithms are being developed to accomodate the size and speed limitations 
of on-board computer systems. One class of these algorithms consists 
of square root sequential filtering methods. The purpose of the square 
root filters is to reduce the likelihood of filter divergence which can 
occur as a consequence of a small computer word length. 

In this initial study, a new method for the time update of the 
square root covariance matrix was developed. In addition, this time 
update method is compared with another square root covariance propaga- 
tion method to determine relative performance characteristics. Compar- 
isons are based on the results of computer simulations of the LANDSAT-D 
satellite processing pseudo range and pseudo range-rate measurements 
from the Phase I GPS. A summary of the comparison results is contained 
in the following paragraphs. 

Summary of Results 

In square root algorithms that employ triangular square root co- 

variance matrices, time propagation by transition matrix methods destroys 

the triangularity of the square root covarian matrix. Retriangulari- 

zation is required after each time update. Such a procedure is required 

T 

of the square root filter (the UDU algorithm) currently proposed for 
the LANDSAT-D computer. In addition to tlie computational burden of the 
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retrlangularization , Lhe effects of process noise can only be approxi- 
mated, if the algorithm is to be computationally efficient. As a part 
of this study, an algorithm has been developed which Integrates the 
square root covariance directly in its triangular form. Retriangulariza- 
tions are not necessary in this algorithm, and the effect of process 
noise is Included exactly. 

Appendix A contains a der'vation of the proposed propagation algor- 
T 

ithra for the UDU filter, as well as the results of a performance compar- 

T 

Ison between the proposed method and a UDU algorithm using a transition 
matrix time update. Two versions of a standard formulation of the Ex- 
tended Kalman Filter (EKF) are included in the comparison, also. The 
results show that, for this test problem, the proposed propagation method 
has superior performance to tlie transition matrix update formulation in 
terms of efficiency and accuracy. Its efficiency is only marginally less 
than that of the EKF formulations. Detailed results are available in 
Appendix A. 

The results of additional algorithm comparisons are contained in 
Appendix B. Two more square root methods, the Potter and Carlson algor- 
ithms have been included in the comparisons. A directly integrated square 
root covariance propagation algorithm, similar to that derived in 
Appendix B show, again, that the direct square root covariance updates 
can be competitive with the transition matrix formulations in terras of 
computation efficiency and estimation accuracy. 
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Nomenclature 

a serai-major axis of sate Hi to 

orbit 

a^ atmospheric drag acceleration 

b filter model clock bias 

bjj^ satellite clock bias 

b GPS clock bias 

s 

c speed of light 

d ballistic coefficient 

c eccentricity 

f true anomaly 

g gravitational accei oia t i on 

h satk-'llite altitude 

li scale height for density nunlel 

o 

i inclination 

k density model scaling factor 

n filter model clock drift 

n^ satellite clock drift 

n GPS clock drift 

s 


r satellite inertial position 

vector 

t true time 

T satellite clock indicated time 

V satellite inertial velocity 

vector 

V , veloeitv vector relative to the 

rel , ' 

atmosphere 

V , magnitude of v , 

rel rel 

V range measurement 

V" i ini’i -rate ir.ca:nn'ir„eMt 

('oi r 'lation parar.eter for clock 

i::' 'tb> ’ 

atiao.sp'aeri c density 

^ atmospheric density at reference 
aJtitudt! 

; :',i 'on'c t r i c rang.a 

.1 geometric rang, e-rate 

M argument c'f jiericenter 

: longitude of ascending node 


Abstracj;. 

A method for propagating the square root of th.e state error covariance matrix 
in lower triangular UDU^' form is described. The projiagat i ('n method can bo 
combined with the UDu'^ measurement inct)r]iorat i on nlgoritlim to obtain a complete 
square root free triangular estimation algor itlim. The method is compared with 
(1) the UDU’^’ state transition matrix propagation algorithm and (2) the con- 
ventional sequential estimation algorithms on the basis of estimation accuracy, 
computational efficiency and storage requiremi-nts, by a simulation of LANDSAT-D 
processing data from the Phase 1 (ilobnl Positioning System. The nuniericni 
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rssuXes iEdlcat® Chut, while slower fhan the ennventiooal asthods, the 
Method proposed here is aore efficient than the previous UB factorlsaclone 
with regard to faBputer storage and comptit.it I <m lime, and leads to the 
most accurate estimate of any of the methods considered. _ 

Introduetioo 

Square root filter formulations have been propo.sed as ;3 means of 
elitfiinating the problem of filter divergence in the real-time application 
of sequential estliuiitlon algorithms. In these methods, the state error 
covariance matrix is replaced by its square root during the propagation 

1 « (C 

and update of the estimate. The state error covnrl.nnce matrix does 
not appear explicitly and, if it is required, it can be obtained by 
multiplying the square root covariance by Its transpose. Consequently, 
it will always be simi-posltive definite. 

In the initial f ormulationsi ** the enhanced numerical stability 
was obtained at tlie expense of lncrea.sed computation complexity, and an 
associated increase in computation time. Tn RoT. 5, a square root 
measurement update method is proposed which offers potential improvement 
in the computational efficiency of the square root filtering methods. 

This efficiency is created fay maintaining the square root covariance 
matrix in triangular form. Following a procedure based on Givens's 

transformation® an algorithm has been pro])oscd^ wltich factors the st.ite 

T 

error covariance P into the form P = UDU , wlicre U is unit upper triangular 

T 

and D is diagonal. Using the UDU fact or i :>.n t i on eliminates the square 
root functions present in the algorithms discussed in Refs. 1-5. The 
measurement incorporation formulation derived for this factorization 
technique is summarized in the Appendix. The proposed algorithm’ for 
propagating the estimate is summarized in the following paragraphs. 


Tlie discrete time propngntlcm oqimtlon I'nr the Ktatc error co- 


variance matrix Jh 


y\ 

where P, , , is the a priori state error covariance matrix at t. , , , P, is 
k+1 k+1 k 

the a posteriori state error covariance matrix at tj^, ^^^k+l’^'k^ 
state transition matrix used to map the state from t^ to and 

is the matrix \^^hich accounts for the effects of process noise in the 
interval from t^^ to The matrices 4>(t^_j_.j ,t^) and satisfy the 

following equations: 


( 1 ) 


4>(t,t^) « A(t)^(t,tj^) ;(f(tj^,t^) = T 


k+1 


k+1 T 

<l>(t^^l,T) Q(T)<i>^(tj^_^^,T)dT 


( 2 ) 


( 3 ) 


whe*-e A(t) is a known nxn time-dependent matrix and Q(t) is tlie process 
noise covariance matrix. 

The discrete square root time propagation algorithm, based on the 
T 7 

UDU transformation, can be summarized as follows : Forai the two 
matrices 


\+l = ^^^‘^k+1 : \+l^ 


( 4 ) 




fD, 

: 0 

k 

0 

• V*' 


( 5 ) 


^ 


where 


k+1 


'k+1 




( 6 ) 


The process noise covariance matrix, Qj^ , is assumefl l:o be ccmstant In tlie 


integration interval At S t, ^ - t, , and , 

k + J- k k ‘r 1 

approximated as 


, as defined In Eq,(3), is 


^k+1 '' ^k+] 


k-H 


(7) 


Then, the updated factors (U, ,,,D, are obtained in upper t r iangular and 

K+1 k+ L 

diagonal forms, respectively, by performing a Modified Weighted Gram- 

Schmidt orthogonal ization on the matrix W, ^ , Vvdiere its columns are 

weighted by the diagonal matrix 

The calculation of $ ,tj_) requires the integration of nxn 

equations in addition to the n-state equations. The determination of 

B, ,, necessitates an nxn quadrature. Tlieret ere , the total number of 
k+1 ‘ 

equations to be integrated is 2(nxn) + n. ihe I’DH* formulation proposod 

in Ref. 7 approximates by an analytical trapezoid-rule integration 

which eliminates the nxn quadrature. The error introduced by this aji- 

proximation can be neglected if the propagation intc'rval (t -t ) is 

k+1 k 

small. The effect of error accumulated o\'cr long, red i (.• t i on intei'vals, 

during loss of tracking or dn.ta drop-outs, must bo considered to ascertain 

the accuracy of this approximation. 

The matrix multiplication, I'U , combined voith the creation of the 

k 

augmented matrix » destroys the triangularity of the square root 

covariance matrix. The application of the modified Gram-Schmidt ortho- 
gonalization procedure is required to rctriangularizc the UD factors 
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at tha timffi each meaaurQment Is pracfssftl . Tlif atidi'd computationnl 
burden of this orthogonalization at each observation point could be 
eliminated if the square root covariance matrix were propagated without 
Che loss of its triangularity. 

In this investigation, a method Is proposed which allows the inte- 
gration of the continuous state-error covariance differential equations 

in square root form. The derivation follows the approach used in Ref. 8, 

T 

but Che results are based on the P = UDU decomposition. The new algorithm 
can be combined with a triangular measurement update algorithm to obtain 
a complete square root estimation algorithm for which square roots are 
avoided. In addition, the effects of state process noise are included 
without approximation. 

The Square-Root Propagation Equation.^ in Triangular Form 

The differential equation for propagating the state error covariance 
matrix can be expressed as 

P(t) = A(t)P(t) + P(t)A'‘(t) + (nO (8) 

where P(t) is the a priori state error covariance matrix, A(t) is the 
nxn linearized dynamics matrix, and Q(t) is the process noise covariance 
matrix. Each of the matrices in Eq.(8) is time-dependent in the general 
case. However, for simplicity, the time dependence will not be noted 
specifically in the following discussion. 

If the following definitions arc used, 

P 2 UDU^ ; Q H Q/2 (9) 

and if the first part of Eq.(9) is differentiated with respect to time 
and substituted into Eq.(8), the results can he rearranged to form 
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(05 + M - QU-'^ - aud)u^ + U(UD'' + f-- ifV - dOV) - o. 

4 \ ]• ' J 

Noting that the first term of (10) is the irnnspose of the second term, 
and making the following definition: 


C(t) H (UD + 


- QU ^ - AUD)U^, 


( 11 ) 


one obtains 


C(t) + C^(t) « 0 

Relation (12) requires that C(t) be eftncr tlio null matrix or, more 
generally, skew symmetric. 

Equation (11) can be aimpHficci by sdoctivoly carrying out tlie 
— T -T 

multiplication of the -QU term by U to yield, after terms are 


rearranged , 


(UD + ^ - AUD)u’^ = Q + C(t) C(t) 


(13) 


Equation (13) defines the differential equ.itions for U and D to the 
degree of uncertainty in C(t). Since tlie unknown matrix C(t) is 
skew symmetric, tliere exist n(n-l)/2 unknown scalar quantities in 
Eq.(13). The problem considered here is one of specifying the 
elements of C(t) so that U is maintained In triangular form during 
the integration of Hq,(13). (Tlie de.rivatLon pursued here assumes that 
U is lower triangular and D is diagonal, although an algorithm for an 
upper triangular U can be obtained ns easily.) The following definitions 
are made to facilitate the solution to the problem posed above. 


nn 

T = AUD ; M E UD + ^ - T 


With these definitions, (13) is expressed as 


(14) 


MU = C = Q + C(t) 


(15) 
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Since U and U in Eq.(13) are lower trlanp.ular, and since from (12), 

C(c) is skew symmetric, several observat Jons ran be made regarding 
Eq.(15). Tlicrc are n(n-l)/2 unkown elei.'enis in C. The produets UI) 
and UD are lower triangular creating n(n-M)/2 unknowns. Therefore, 
the nxn system of equations (13) has [n(n-l)/2 + n(n+l)/2] =• n ;c n un- 
knowns which can be determined uniquely. 

An expansion of Eq.(I5) into matrix c-lemonts indicates tlie method 
of solution. 


1— 1 



T j 2 • • • 

-T, ^ 
1 n 

> 

1 

^21 

••• ^nl 
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M 

"22 • • 

-To 

2n 


• 

1 
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\2 

• 

M 
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0 
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21 
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nl 
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^22 • ‘ 
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. 
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; 
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‘'n2 • • 

• h 

nn 


In Eq.(16), Q i.s a.ssumed to be a diagonal mati ix witli elenieiUs 
q^^ " ^11"^^ ’ (Thi.s assuniUioii ran bo gi-no’ a 1 i r.ed 


(16) 




Vi- 


to allow other non-zero terms in the Q matrix with only a slight in- 
crease in algebraic complexity.) Each row of the upper triangular 
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portion of th« C matrix In Eq.(16) 1« detcTitihied as the product of the 

-T 

corresponding row of the M matrix with the appr.- 'r i 'i <• cnlumn of the U 
matrix. After an upper crlangul.ir r!>w of C [,s tiei ^ , n s rei; j (he coiulitum 


from Eq.(12) that C 


- C.. (1-1,. ..,n ; i“l i-l) is invoked to 

i.1 .1 i 


evaluate the corresponding lower triangular eolumn of C. Then a column 
of the lower triangular elements of M can be evalu.ited. Once the 
elements of the M matrix are determined, the next row of the upper tri- 
angular C elements can be computed along with a column of the 0 and 5 
elements. This process is repeated until all U and D values are 
determined. The implementation of this approach proceed.^ as follows. 
From Eqs. (13) and (14) one can write 


M + T - UU + 


(17) 


The expansion of (17) in summation notation gives 


M + T. , * 51 U,, d , + ). — 5 — ^ 

ij iJ ,,.1 Jk kj 1 ^., 2 


i-^1 , . . . ,n ; j = l, . . . , i 


(18) 


But, since D is diagonal, (18) becomes 


M.. +T,. = U..d.. + 
ij ij ij JJ 


U. ,d. . 
-LL-U. 


i=l, . . . ,n ; j=l, , , . ,i 


(19) 


For i“j , U = I and U.. = 0. Therefore, (19) becomes 
, ij ij 


d. . = 2(M. . + T. J i=l , . . . ,n 
11 11 il 


( 20 ) 


For i ^ j« (19) i® rearrnnju’d to obtain tlif d il t or^nt ial ec,uation 


U, ,d 


°ij ■ '"n 


1-2 , . . . ,n ; .1 = 1 i-1 


( 21 ) 


LquaCions (20) and (21) ;irt* (lie i(jrr!s o; tiu' d i 1 ! v-rant ial uquati.ms 
to be employed in the derivative routine oi ,i numerical integrator. Tlie 


elements and are computed as defined in F.q.(14). The pertinent 


equations can be combined to obtain the following algorithm. 


Triangular Square-Root Propa Ration Algorit hm 

Given the elements of the square root slate erri'ir covariance in 
lower triangular UD form, Q 2 Q/2, and A(t) , the differential 


equations and d^^^ can be computed as follows 


T 


. , = A . , U, . d . . 

k=j -11 


i^'-l . ■ . • ,n ; ,'p I 1 


( 22 ) 


C . . = T. u . , - ): T V 

k=l k=i-U 1’^ 


>■ 1 n ; 1= ■: fl , . . . ,n (23) 


i-1 


“ii - ‘'ll ''idik -■ 

k-1 


(24) 


= -C.i - V M.kii,^ i = 2,...,n ; j = l,...,i-l 


(25) 


d. . 
11 


2(M,, + T. .) 
ii 11 


i=l , . . . ,n 


(26) 


ij 


U. .d. . 

(M + T . - -i^-^)/d. . 

ij Ij 2 i-i 


.li 


. ,n 


.1 = 1 ,i-l (27) 





The propagation algoritliin sunun.irizc-d In ( 22 ) Uitnup,h (27) ran In- 

combined with the algorithm for incorpornt I nr. on nb.si‘rvatlon to obtain 
a complete sequential estimation algorithm In which, tin; covariance 
matrices P and P are replaced by the factt'r.s (U,D) and (0,D), respect ' vely . 
The algorithm given in tlie Appendix nstnn.ies t:iar only a .sinp.le sc.ilar 
observation is processed at inndi (diservat ion .po.h; howi\i i , (lie . 1 1 p.or i t Inn 
is applicable to tlie case of multiple ohsi.' r\’.it io:n- , 1 ; riven ■.■noi-h, if 
the observatio' I'rror.s are assiiinod to he uiuor n 1 a t ad . 

Numerical Comparison 

In tiie folli'wing luimericn.l example, 1 !ir i v. 1 n.otlioPs for I'la na.M! ion 
T 

of the UDU factored covariance matl ix arc- cor-.p.iroc; 10 dolermin.- the 
relative computation speed and estimation actniracy. As a basis for 
determining their .absolute performance, nuiroric.il results arc obtained 
with the conventional Extended Kalinan-IUicy filter using both I’(!;;.(l) and 
(^) as the bases for prop.agating the st.ate eriao- .•(' .-.ar iam-e matrix. The 
numerical comparison.s are made by using e.u-h of the algorithms to process 
a set of simulated (llobal Positioning .System ((IPS) r.inge .and rang,e-r.ate 
observations obtained by tiie 1,AN1)S;\T-D sp.'ivcr.ii't . A det.iilod discu.s.sion 
of the GPS and the associated n.ivig.ation mo.a.sr I'ericnts is i'.iven in Pei. i2. 
Since the range measurement will require a precise r.ie.isureme,';t of the 
time interval between ru'gnnl transmission ' l om o-ne of the (P’S .a.a t e 1 1 i t c.s 
to reception at the I.ANUSAT-P sp.aoeer.afl , ( T,< cl..o!, oitoi- ran.': : 

modeled and included ns part of tin- over.:l; .‘-tato vecLor, !'hu ])riinar;.' 


clock errors are the bias and drift. 
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The dynamic model used lor 1 lie motion ol I.ANDSAT and Llie associated 
model of the satellite's clock are coirbinod to obtain a filter model which 
contains nine state variables. The nine cimiponents of the state vector 
are: position (r; 3x1), velocity (v; 3x1). dock bins (b) , clock drift 

(n) , and clock-drift model correlation pa rane t ir (P). The differential 
equations defining these parameters in the filter are: 


r *= 

V 

(28) 

• 

V >» 

S + + ^v 

(29) 

b - 

n 

(30) 

• 

n “ 

Pn + r, 

n 

(31) 


' P 

(32) 


The stochastic processes, ^ , C «^nd f, arc white noise forcing 

v n 

functions which are assumed to have the ftil lowing statistics: 

E[C(t)]^ - 0 ; E[f,(t)f;'^'(T)]_^ = (!.U)S(t-0 (33) 

where tiie subscript i indicates the apnropr in te member of the set {v,n,Pf 
and 5(t-x) is the Dirac delta function. 

For computational efficiency, the filter assumes -.Imple nujdels for 
the gravitational acceleration, g, and fnr ibe a i mospher 1 c drag .icce U-rat i on , 
a^. The geopotentiaJ model adopted for the fill.er is obtained l)v trun- 
cating the Goddard Eartli Model (GEM7)*' to the lourth degree and ordiU'. 
d’O drag acceleration is calculated as 
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clYv ,v , 
rol ri'I 


(3-'0 


where tliu atmospheric liensliy, y , Ik approximated by the exponential 
model : 


Y 


I 


-k(h-h ) 

V a 

o 


i '•)) 


The values of the dr.i)\ moilel pai .miet ei:; a; aii",.,! ioi tliis iiivesi i I'.ai 1 on 
are! h •• tUiO.OOt) in, y - ''./a x 1(1 ^ ' kr./n; , k '.SS x 10 m, aiul 

O i' 

d 1.18 X 10 ‘ m ’ /kp. . 

The estimates of the clock bias, b, aiul arift, n, arc used tv) 
pri'dict the true time of the user's clock, t, bv the equation 

t - T - b - n tt-t ) (16) 

o o o 

where the subscript (o) imlic.iti-s the epoch <’i t !;c last estimate of 

the p.'iramet or s , The (luantity T is the t inu' .is indicated by the LANDSAT-D 

c 1 o c k . 

The observations used ior this study arc p.soudo ranpc and psciulo 
r.tnpe-rate measurements as ('bsorved by the 1 AN'IiSAT-P satellite usinp 


the si.x satv’llites ol 

t he 

Phase 

1 OP.S ci'iis 

I e 1 1 a I i on . ■ 

Tbe computer 

s<0'tware used to s 1 mu 1 

at i- 

1 he ob 

servat i oin; 

as wo ! 1 as 

I ni l her detai 1 s 

on the simulation proc 

■edur 

0 ail' 

d i scussed 

in Kc 1 , 1 (' . 



The models used to I'.cncr.ito the siiailjyc.l ri a: nia men! s. ha\-e the 

t oi m : 

Y ■ i (b,, - b ic I /. (-17) 

P Is p 

Y' . I' I (b, - u )c I • (’38) 

p ^ S l' 

vlu^vo p and p arn Iho l nu’ vaiuns i .iiu’ r.in/,r* i .itn lu'LWLU'n 

1 .ANnS A*l “O aiui a p,i\’nn td'S s.i i t* I 1 I I t* , V .ni.i V* .irt- tin nu'.iaurt'tl 


values of the ratn*e and range-rate; bj , ii(, , b^, 

and drifts In tite salelllli' i-lnrk ami in t!i<' (b'S 

1 ? 

and c is tlie speeci of ligliL. 

The measurement Y anti Y* is proei'sseci i^v 
P P 

filter using ti\e model 


and n are tlie Ijiases 
s 

e locks, respec’ L i ve 1 y ; 


tlie I.AND.'iA’I'-D navigation 


Y = p + cb + C ( 

P P 

Y* = p + on 4- r, ' (' 

P P 

That is, in the filter model the GPS clocks are assumed to be perfect 
and the total time error is assumc'd to lie contained in tlie I.AMDSAT-l) 
satellite clock. 'I'he interval between observations 's assumed to be 
six seconds and observations from only one G]’;: satellite can he obt.-timui 
during any six-second interval. The. obse.rv.at i ons from the visible 
satellites are processed sequentially. 

Tlie GPS satellites are assumed to be in circular orbits (<.'-=0) about 
a point mass earth with inclinations of 6J" and periods of 12 hours 
(A3, 200 sec). Three satellites arc equally spaced on each of two 
orbital planes. The orbital elements are rt- f eraiiced to a coordinate 
system whose, xy-plane i.s the earth's equator and whose ;<u-plane lies 
along the Greenwich meridian. 

The epoch condition fur LAMDS.AT-D was chosen so that the resulting 
simulated observations would nccurntoly reflect the possible extremes 
of GPS satellite visibilit.v. Tlie epoch, Ciem,,,^^ chcsi'n .ire 
a E 7.086901 x 10^’m, e H .0001, 1 9S:i8!, . .3TV:878, oj 180?, 

and f(true anomaly) 2 -185°. Thu elements aiu' specified at a GPS svc.te 


time t=C. The. epocdi eli'iiiicnts for GPS arc ,s; 


ec i 1 ii d at a sveitcm time of 


-7200 sec. The difference In initial c]>o<h s is included in the filter 
program's update of user and GPS states. 

The values of tlic LANDSAT-D and GPS clink pliase and ircquency 
errors are simulated .as the sum iif tlirec d liferent error sources; 
a noise-free phase error with a polynomial form, an error due to ex- 
ponentially correlated frequency noise, and a random walk bias error. 
Tlie exact form of the error models and tlie coot I i e ient s used in the 
models are given in Ref. 10. 

The numerical simulations of the filler performance were made 
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with the following initial conditions: 


The off-diagonal terms of the state-error covariance and noise covariance 
matrices are set to zero initiallv. 
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The data in Table 1 and Table 2 show ilu- relative performance of 
the four algorithms when each processed a simulated ohsiTvat inn data 
set with a duration of approximately lODOD sec. The five columns in each of 
the tables are, respectively: (1) the tot.il (H’l’ lime required to perform 
measurement updates of the state and covar i aiu’e , (2) the total CPU time 
required for propagating the state and eovar i aiue , (3) the total CPU propa- 
gation time normalised by the fastest CPU proiuigat iiui time, (A) the total 
CPU time required for proptigation and measurement updates (the sum of 
columns (1) and (2)), and (5) the RMS of the position magnitude errors 
and velocity magnitude errors over the duration of the simulation. All 
CPU times are listed in milliseconds. Position errors are given in meters 
and velocity errors are given in motors per second. The algorithms are 
ranked in the tables in order of increasing total computation time. 

Time propagation is ])erformed with a fixed-step modified F.uler 
Integrator, which requires two function eva 1 u.i t i on per step. The. inte- 
gration step size for the data of Table 1 is six seeonds, equal to the 
time interval between CPS observations, I'lie data of Taiilo 2 i osu ! l from 
integration with a three-second sti'p si.xe. 

'J’he relations for the propagation of Llie time bia.s and drift as 

approximated by Eqs. (30), (31) and (32) liave simple analytic sc'Iutions. 

This allows certain elements of the stat.' transition matrix tc' be updated 

analytically. The implementation of the riH'(v'> ■uu! fiKC (i) .ilgoritlims 

has taken advanttige of this slmpl if It-at ion to reduce the number of numerically 

•) 

integrated differential equations below t lie t lieori't i ca 1 v.aUu' of n" . I'he 
high degree of coupling in the covariance diflerential equations for the 
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(P) und (U,D) algorichma does not permit a convenient reduction in the 
integration vector aiaa. A total of n(n+l)/2 covariance equationa haa 
been integrated in the siwulations descrlbc'd hern. 

Wliile the programming effort required to implement the (U,D) formu- 
lation will be greater because of the recursive nature of Eqa. (22) through 
(27), fewer computer storage locations will bo required to execute this 
algorithm since the total number of equations involved in (22) through (27) 
is n(n+l)/2. This value compares with the (nxn) computer memory locations 
required to Integrate and store the equations. Siiice there are some 
zeros in the | equation, one can reduce the stor.nge requirements of this 
method at the expense of added programming eoinplexity . Further comparisons 
of the computer storage location requirem .its for the.se two algorithm.s are 
given in Ref. 10. 

The numerical results shown in Table 1 indicate that the (U,D) algo- 

X • 

rithm is competitive with botli the UDU (>J>) and conventional formulations 

in terms of CPU times and estimation accuracy for tliis filtering problem. 

llie (U,D) method is faster than the UDU (t) algorithm for the six-second 

integration interval. Its position estimation error is lower tlian that 

for any of the other algorithms. For the three-second .step size results 

in Table 2, the (U,D) algorithm remains competitive in terms of estimation 

T • 

accuracy, but is no longer as fast as the UDU (!') algorithm. With the; 

decrease in integration step size, the number <'f expensive (U,D) function 

evaluations has increased, lujt tlie luiml'er oi’ lime consuming or tliogona 1 i - 

'i' • 

zations in the UDU (^^) algoritlim remains the same. This factor causes 

T • 

the UDU ($) algorithm to have a faster computation time, In this case. 
Again, for Che three-second results, the (U,D) formulation yieliis the most 
accurate position estimate. 
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Both UDU methods require higher total computation times than the 
EKF algorithms. This is not an unexpected result as squ.ire root filter 
formulations often incur computation time penalties as the price for 
increased numerical stability. 

These numerical results were generated on a CDC6600 computer system. 

The relative performance of the algorithms will vary ns a function of the 
computer system being used, the dynamic model assumed by the filter, ch® 
method and order of numerical integration, and the integration step size. 

The influence on performance of the latter two factors can be seen in the 
numerical results given in Ref. 10. The (U,D) method becomes more efficient 
as the total number of function evaluations within a given integration 
Interval is decreased. Therefore, the choice of the time propagation 
method should depend on the formulation of the* specific problem under 
consideration and on the supporting algorl ttims and computer system used 
to perform the calculations. 

Conclusions 

Based on tlie results presented in the previous discussion, it is 
concluded that, for the example problem considered hero, the (U,D) algorithm 
is more efficient than the (U,D) algorithm based on the ^ propagation. 
Furthermore, the estimate obtained with the (U,fi) formulation was more 
accurate tlian the estimate obtained l)y eitlu’r the conventional F.KF esti- 
mation algorithms or the (U,D)-^ algorithm. The performance of the algorithms 
will be dependent on the computer architecture and software, the dynamic 
model assumed for the filter and the method used to perform the numerical 
integrations, and will vary as these factors change. 



The measurement update algorithm for the IDU factorization has 


the following form. Using the obsc-rvation Y, , , = G(X, . , , t, 
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Compute residual: 
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Calculate gain and update state: 
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Limited word size in contemporary microprocessors 
causes numerical problems in autonomous satellite navi- 
gation applications. Numerical error introduced in 
navigation computations performed on small wordlength 
machines can cause divergence of sequential estimation 
algorithms. To insure filter reliability, square root 
algorithms have been adopted in many applications. 

The optimal navigation algorithm requires a careful 
match of the estimation algorithm, dynamic model, and 
numerical integrator. In this Investigation, different 
representations of these elements are evaluated to de- 
termine their relative performance for satellite navi- 
gation applications. Numerical simulations are conducted 
using the Phase I GPS constellation to determine the 
orbit of the L/\NDSAT-D satellite. Numerical comparisons 
are made of various square root filter formulations, and 
their dependence on the order of the integrator is 
examined . 


Nomenclature 


a, atmospheric drag acceleration 

d 

b satellite clock bias 

c speed of light 

d ballistic coefficient 

g gravitational acceleration 

G computed range measurement 

Gp computed range-rate meas . 

h satellite altitude 

h scale height for density model 

o 

k density model scaling factor 


n satellite clock .irift 

r iin.Ttial ’,'oalt;ion vector 

t true time 

T satellite clock indicated time 

V inertial velocity vector 

V , magnitude of 

rel rel 

V , velocity vector relative to 
rel 

the atmosphere 

correlation parameter for 
clock drift model 
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The 


Y 

variable atmospheric density 

• 

P 

c 

computed geometric range-rate 


atmospheric density at refer- 
ence altitude 

Pm 

• 

m 

^Pm 

range measurement 
range-rate measurement 

P 

geometric range 

range bias due to error in 

P 

geometric range-rate 

user satellite and GPS clocks 

Pc 

computed geometric range 

Ap 

m 

range-rate bias due to error 
in satellite and GPS clock 


1. Introduction 

The use of artificial earth satellites for accurate dissemination of 
time and frequer y holds high potential for the development of an accur- 
ate and reliable autonomous navigation system [1]. Current satellite 
systems have demonstrated the ability to obtain global positioning of 
points, fixed on the surface of the earth, to an accuracy of two meters 
in three dimensions [2]. While the position error achieved by a dynamic 
navigator, i.e., one moving with respect to the earth's surface, would be 
considerably greater than two meters, this investigation indicates the 
potential inherent in satellite navigation methods. 

The Global Pcsicioning System (GPS) [3], which is being deployed currently, 
is designed to allow a user to satisfy real-time navigation requirements 
by the calculation of position and velocity using simultaneous pseudo 
range and pseudo range-rate measurements from several GPS satellites [4, 
5,6]. The requirements for determining the orbits of low altitude satel- 
lites in near real-time, coupled with the need for increased accuracy, 
generates an interest in evaluating the GPS as a means for satisfying 
satellite orbit determination requirements. With the development of com- 
pact low-power computers and atomic clocks, the ability to perform the 
satellite navigation function on-board the spacecraft in an autonomous 
navigation mode is an attractive alternative to telemetering the GPS 
range and range-rate measurements to the ground for processing by a 
ground-based orbit determination program [6]. 

Allowable computer storage and execution times will place constraints on 
the model and the algorithms which can be selected to estimate the satel- 
lite's state. To minimize the storage requirevients and achieve a real- 
time state estimate, the estimate of the satellite orbit will be performed 
on-board, sequentially, using a Kalman-Bucy filter [7]. One problem 
which must be considered if a sequential data processing method is used 
is the problem of filter divergence [8]. The divergence occurs due to 
either (1) dynamic or measurement model error, or (2) numerical errors 
introduced during the computation process. Since most computers for 
autonomous satellite navigation will have a short wordlength, this second 
cause of divergence will be of considerable importance. 

The problem of filter divergence has led to a number of s.'.udies aimed at 
the development of stable estimation algorltlims . The square-root measure- 
ment update algorithm proposed by Potter [9,10] has been used in a number 
of applications to prevent divergence caused by a computed non-positive 
definite covariance matrix. The algorithm proposed by Carlson [11] 
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allows ihc rooasurcmotir up>iato to bo .iv-oomp I i :!ioil in .in otfioiont nunnor 
by nuiintaining a triangular square root lovarianoo r.iatrix. In (!2), a 
triangular decomposition which ia free of square root operations is pro- 
posed. The "square root free" algorithm is reported Co be taster chan 
the Carlson-Cholesky algorithm. The squ.ire root covariance matrix must 
be maintained in triangular form during the time propag.itiun it these two 
algorithms are to be .ippllcable. Usuallv, the squ.nre root co'-.'r iance 
matrix time update destroys the matrix t r i.inc.u l.u i ty ami .i ret i iangu lar i- 
xation at each potential observation epoch ;;ui.sc h>' emploved. This pro- 
cedure requires an .ulditlonal computation effort iml an .is.soe laced 
computation time pennlcy. A recent devolopmonr bv Tapley, ut al [13, 1'^], 
in which, the time update of the square root oovari.mce matrix i.s main- 
tained in lower triangular form throughout llu> entire estinatlon process, 
can be adopted to obtain a complete tri.mp.ular squari^ root estimation 
algorithm. 

The objective of this i nve.s t i gat ion is to ev.iiuatc che perform.ince of 
the various square root algorithms in performing on-board satellite 
orbit determination using Che Global Posi t ioning .'■'ystera. The evalu.i- 
cion is based on a comparison of che comput.ition time and the .state 
estimate accuracy. The eftect of the numerii'.il inf egr.ic ion method on 
Che estimate accur-icy is con-sidered also. 

2. Filter Model 

Tlie moviels of the .satellite dvnamlcs .ir.l of the observ.'U ion st.iti; rela - 
tion hvive a cricic,al Impact on the compviter stor.igo requirement and on 
the execution time for the n.fvig.jt ioi\ .ilgori i !im. If the models .ire too 
complex, unaccepcablv large storage requi renont s ,uul vU'mput.ic ion times 
will occur. However, If the modi-l.s .iro r-'o simnlo, ch,' .icoui'aov of the 
navigation estinute will lu' dogr.fdcd boiow ar, .i.copt.ib 1 o v.iluo. !'hc 
following model Is selcctcvi .is ,i v-omp ror. i sr botvocu the roo,u i ivmont s fur 
accuraev and efficiency. 


Dynam ic K quati ons. Tlie dyiKuni,- model f.-r '.ho moti.-u of ; !\o u,-;er s.itelli.te 
and the associated model of the satellite .-look behavior .ire eombined to 
obtain .a filter model which contains nine ate v.iri.iblcs. Tlie nir.e com- 
pov\ev\C.>5 of the state veeior are: posit i.ai U: ' ' P, vel.n-ity vv; > ■- 1 ) , 

clock bias (b), clock drif t fiO, .nui cl,'ok-dri: t ■.node’ c,irre,l.it ion j).,ira- 
meter (3). The viifforentl.il equ.'ition.s vietining these p.ir.imeters in the 
filter are: 

r ■ V f 1) 


g •+ .1 , + t, 
d V 


pn 


n •• Bn + K 

n 


i, -j ) 


« 

where ( ) •• d t. 
clock time. 


),.'ilT, c.g., the ;ndepo:ivii nl 'ca 


fS) 


able i.s t!\e satellite 









The stochastic processes, and arc white noise forcing functions 

which are assumed to have the following statistics; 


E[^(t)] - 0 

ElCCOc^CT)]^ . Qj^(t)6(t-r) 


( 6 ) 


where the subscript i indicates the appropriate mernber of the set 

{v,n,B). 


The filter uses relatively simple models for gravitational force, g, and 
the atmospheric drag forces, a^. The geopotential model adopted for the 
filter is obtained by truncating the Goddard Earth Model (GEM) 7 [15] to 
the fourth degree and order. The drag acceleration is calculated as 




(7) 


where the atmospheric density, y , is approximated by the exponer/.ial 
model; . /u u \ 

Y . Y (8) 

o 


The values of the drag model parameters assumed for this investigation 
are: h^ •• 840,000 ra, Y “ 5.74xlo~^‘‘ kg/m^, k *■ 7.58'<10“®m, and 

d - 1.18 x 10~^ mVkg. ° 


The clock bias, b, and drift, n, estimated 'iy the filter, are used to 
predict the true time of tl'.e user's cloc’.. bv th.o equation 

c “ T “ b - n ( t - t ) (. V ) 

o o o 

where the subscript (o) indicates the ep.'cb. >; tl'.e last ova ! uat i (,'u of the 
parameters . 


Varlationa.l Ec uat ions ■ Linear variation.'.! oer.ati^ns are required to 
implement the sequential estimation algorithm [8]. The equations are 
derived by the linearisation of Eqs. (1) through (5). Furtltor details 
can be found in [13]. 


Observation State Relation. The observation type.s processed during the 
study are pseudo range and pseudo range-rate. The actual measurements 
can be expressed, mathematically, as: 
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(10) 

(11) 


where p and p are the measured values of tlie range and range-rate, p 
• « 
and p are the actual geometric range and range-rate, Ap and Ap are range 

m m 

and range-rate biases due to errors in the user's clock and the Gl’S satel- 
lite's clock, and the are gero moan white noise sequences with knot/n 
variances. ^ 
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The modeled values of the measurements have tlie f*7rm 


G =p +(b~b)c 
pc s 

G* = p + (n ” n )c 
pc s 


( 12 ) 

(13) 


where p^ and p^ are Che computed values of range and range-rate, b, n, 

b and n are tlie predicted biases and diilts in the satellite clock 
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and in the GPS clocks, based on previous moasuroTnenCs , and c is the 
speed of light. 


The linearized observation-state matrix H is computed by caking partial 
derivatives of G^ and with respect tc the state. The complete ex- 
pressions for G and G*, along with tli<> partial derivatives, are .;i.v(>n 
in [15]. 


3, Observation Simulation Model 

Observations generated for this study arc pseudo range and pseudo range- 
rate measurements as observed by the LANDS.-\T-D satellite using the six 
satellites of the Phase I GPS constellation. Ihe computer software used 
to generate the observations is a modification of the program developed 
by Kruczynski [5]. Further details on Che software modifications, which 
were made to simulate the orbit of Che L-\XDSAT satellite, are given in 
[15]. 

llie simulation philosophy was Co produce a physically realizable sec of 
data points against which the filters could be tested and evaluated. 
Simplified models of the GPS satellites were used to reduce computer time 
requirements . 

Tlie description of the observation generation program can be broken into 
four basic areas: (1) simulation of GPS satellito motion, (2) simulation 

of LANDSAT-D motion, (3) simulation of clocks, and (4) simulation of the 
measurement process. 

Simulation of GPS Satellite Moticn. For s imp I ic i tv , tl.c GFS satellites 
are assumed to move in circular orbits around a p(Unt m.tss eartli. T'ac 
GPS satellite motion can be determined, tlicn, by using a closed form 
solution and a set of Keplerian elements defined at some specified epoch. 
This approximation will not have a significant impact on the results pre- 
sented here because of the relatively short time interval involved in 
the simulation. 


The epoch orbital elements for the Phase 
following : 
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The GPS satellites are assviaed to ha la circular orbits (e * 0) with 
iacliaatioas of 63* and periods of 12 hours (43,200 sec). Tliese ele^ats 
are referenced to a coordinate system whose Ky-plan® is the earth's 
equator and whose Ke-plane lies along the Greenwich Meridian. 

Simulation of LAWDSAT~D Motion. The force model used to simulate the 
LANDSAT-D motion contains the effects of the earth's non-spherical mass 
distribution and the effects of atmospheric drag. Tlie gravitational 
accelerations for the observation simulations are obtained using the GEM7 
geopotential modal truncated at order and degree 8. The drag acceleration 
is computed by using Eqs. (7) and (8), with the same set of constants as 
specified previously. 

The epoch conditions for LANDSAX~D and the GPS satellites were chosen so 
that the resulting simulated observations would accurately reflect the 
possible extremes of GPS satellite visibility. The epoch elements chosen 
for LANDSAT~D are; 


a 5 7.086901 x 10®m 
e = 0.001 
i = 98?181 
n = 354?878 
ui = 180“ 

f (true anomaly) = -185? 

These elements are specified at the GPS system time t =■ 0 , The epoch 
elements previously listed for the GPS constellation are specified at 
the GPS system time t = -7200 sec. This difference in initial epoch is 
included in the simulation program's updates for the user and GPS states. 

Simulation of Clocks. The values of the L/INDSAT-D and GPS clock phase 
and frequency errors can be calculated as the sum of three different 
error sources: a noise-free phase error with a polynomial form (e^) ; 

an error due to exponentially correlated frequency noise (e^) j a 
random walk bias error (e^) . 

The polynomial error term is expressed mathematically as: 

£^(t) = a, + a^Ct-a^) + (a^/2) ( t-a^^) ^ (14) 

The exact form of the error models and the coefficients used in the models 
are given in [5]. 
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ilaulftElcm of Msaaurasenc Procsas. To generate the obtarvationi , the 
CPS satellites are "sampled" every six seconds in ascending numarical 
order. Upon calculation of the geometric range vector between the 
LANDSAT-D end a GPS satellite, a measurement is accepted if the GPS 
satellite is no more than 20® below Che LANDSAT-D local horizontal. If 
a satellite is rejected because of this geometric constraint, the next 
higher numbered satellite is si»plad in the same manner. The procedure 
is repeated until a "visible" GPS satellite is found. If non® of the 
six satellites is considered visible at a particular time, no measurement 
is taken. Then, all satellite vehicles are propagated forward six seconds 
and the procedure is repeated. A random number of measurements is re- 
jected in an effort to simulate measurement losses due to actual system 
problems such as periodic failures in signal acquisition or bad data in 
the GPS transmission. A Gaussian error term is added to each of the 
geometric observations to account for purely random anomalies in the 
measurement process. Tlie standard deviations for the random errors are: 

Op ■ 2. m, and * .2 m/sec. If no satellites are visible at a particu- 
lar time, the user's clock errors and position and velocity magnitudes 
are recorded on a file to be used to compute navigation errors during 
data gaps. 

4. Filter Algorithms 

The filter used to process the GPS measurements will consist of two major 
segments. These segments are: (1) the measurement update segment, and 

(2) the time propagation segment. The measurement update segment receives 
the observations at a given time epoch and processes these observations 
to obtain an updated estimate of the state. The propagation segment maps 
the estimate and the associated state error covariance matrix forward in 
time to the next observation epoch. For each measurement update algorithm, 
there are two propagation algorithms which can bo considered. Tlie primary 
difference in the propagation algorithms is determined by whether one 
integrates the state transition matrix for propagating the state error 
covariance matrix or whether the differential equation for the state er- 
ror covariance matrix is integrated directly. The filter algorithms 
compared in this investigation are: 

(1) The Extended Kalman-Bucy Filter [8,10] 

(2) The Carlson-Cholesky Filter [11,13] 

(3) Hie Potter Filter [9,10] 

(4) 'n-ie UDU Filter [12,14] 

The state error covariance matrix can be based on the following sec of 
differential equations: 

P(C) * A(t)P(c) + P(c)a’^(c) + Q(t) (15) 

where A(t) = [3F(X,t)/3X] , the nxl state vector, X(t) is defined to 
have Che components X (c) = [r '.v !b,n,;-], F(X,t) is an nxl vector whose 

ik 

components are the right-hand sides of Eqs. (1) through (5), and [ ] 
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indicates Chat the elements of the matrix are evaluated on the reference 
solution X*(t), Alternately, the covariance matrix can be propagated by 
using Che state-transition matrix, t(c,tj^), where (JCtjCj^) » A(t)$(t , Cj^) ; 

<I>(tj^,tj^) “ I. ITic integral of Eq. (15) can be expressed, then, as: 


\+l “ ^ ^'^k+l*'^k^ 


■k+l 


'J>(t,T)Q(x)<I' (t,T)dT (16) 


Rather Chan evaluate the integral in Eq . (16), an average value, , is 
used, where 


- avg 


'k+l 


<I>(t,T)Q(t)‘l> (t,T)dT 


Using Eq. (17), Eq. (16) can be expressed as 


^k+1 * '^^‘^k+l’‘^k^^k ^‘^k+l’‘^k^ ^k 


(17) 


(18) 


The square root estimation algorithms can be based on either Eq. (15) or 
Eq. (18). In the numerical simulations described in the next section, 
methods based on both approaches are compared. The use of Eq. (15) al- 
lows a triangular square root factorization for the covariance to be 
maintained during the propagation interval. Propagation with Eq . (18) 
will destroy triangularity and, after the propagation interval, special 
computation techniques are required to obtain a new square root covari- 
ance matrix. For the square root propagation algorithms based on 
Eq. (15), the relative advantage of maintaining the covariance matrix 
in triangular form is offset by a more comp licated form for the govern- 
ing differential equation. 

Further discussion of the algorithms as well as the specific implementa- 
tion used for this investigation is given in [15]. 


5. Algorithm Comparison 

The numerical performance of the algorithms discussed in Section ^ was 
compared by conducting a series of computer simulations in which the algo- 
rithms were used to process GPS range and range-rate obser\’ations . Ob- 
servations were simulated for a GPS system time interval of 10,000 sec. 
from the LANDSAT-D orbit. This is approximately 1.7 revolutions of 
LANDSAT-D. A history of the number of GPS satellites visible from LAND- 
SAT-D, versus time, is shown in Fig. 1. It can be seen that the number 
of satellites visible varies from zero to six, the number of satellites 
in the Phase I GPS constellation. 

Tlie numerical results obtained with each filter are very similar to those 
shown in Figs. 2, 3, and 4. The plots in these three figures are, 
respectively, RSS position error versus time, RSS velocity error versus 
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tirae, and the error in the ti®e estimate versus time. This solution was 
generated with the UDU(*p) algorithm and a modified Euler Integrator with 
an integration step size of six seconds. 

Two types of error growth are visible in each of the! figures. First, 
there are two tirae periods in each of the plots during which large error 
growth occurs. Comparison of the three figures with Pig. 1 shows a co- 
incidence between the occurrence of the large error and the periods when 
fewer than four CfS satellites are visible. The large errors result 
from the inability of the filter model to predict the state estimate 
through time periods of limited observation data. 

The second type of error growth occurs gradually. After each period of 
low satellite visibility, the navigation error is reduced to a lower level. 
However, the average value of the lower-level navigation error increases 
during each successive period of good satellite visibility. Additional 
simulations run with algorithms other than the one used to generate Figs. 

2, 3, and 4 produced nearly identical error plots. It is likely that the 
long-term error growth is caused by the influence of geopotential model 
error on the estimate of the LANDSAT-D clock parameters. The large errors 
result from the periodic reduction of GPS .-jatellite visibility. Both er- 
ror types require further study; however, since the primary purpose of 
this investigation is the comparison of the algorithm efficiencies, in- 
depth study of the causes of the large error growth and methods for their 
removal is deferred to a later investigation. 

Operation Count Comparison. As an evaluation of the theoretical compu- 
tational efficiency, operation counts of the number of numeric opera- 
tions have been made for the seven different time update algorithms 
under consideration. The operations recorded are the additions (sub- 
tractions), multiplications, divisions, and square roots required to 
perform a single time update of the state error covariance matrix or 
the square root covariance matrix. The operation counts required to in- 
corporate the measurements are discussed in [12]. 

Table 1 and 2 give the number of operations for the covariance time 
update of a system with an n-dimension state vector, whose complete nxn- 
transition matrix is obtained by numerical integration. The dimension 
of the system's state noise covariance matrix Is m. The counts are 
broken into three groups depending on whether the operations occur once 
per time update interval, once per integration step, or once per inte- 
gration function evaluation. 

The operation count per integration step depends on the type of numeri- 
cal integrator being employed. Tables 1 and 2 give the values for a 
second-order Euler integrator with two function evaluations per step. 

The coefficients of the n^ terms in the table will increase significantly 
for higher order integrators (roughly as a factor of k(k+l)/2, where k is 
the number of function evaluations). 

If the assumption of an Euler integrator is maintained and the integra- 
tion step is specified, then the three count groups can be combined to 
give the total number of operations for a covariance tirae update. 
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Aasuming that a time updats occurs avary six seconds and the integration 
step size is six seconds, the total number of operation counts can be 
computed. The result is shown in Table 3. 

The total counts show that direct covariance integration algorithms have 
fewer operations if all transition matrix elements are numerically 
integrated in the (#) algorithms, In this situation, the direct inte- 
gration methods have fewer equations to numerically iategrate, fewer 
operations in each function evaluation, and are not forced to retriaagu- 
larise at the end of the update, 

However, it is often the case that the solution to soa® of the eleffients 
of the transition matrix can be obtained analytically. In such a case, 
the number of terms to be numerically integrated can be reduced. A 
similar reduction in integration vector size has not proven feasible for 
the algorithm based on direct integration of the covariance matrix. 
Therefore, a full n(n+l)/2 set of covariance elements must be integrated 
numerically . 

The actual numerical operation counts for each time update algorithm are 
shown in Table 4. These values reflect the following set of assumptions; 

(1) The state vector size, n, is 11; the measurement noise co- 
variance vector, m, has the dimension of 8. The eleven 
state filter is the maximum filter size used in these studies. 
For the 11-state filter, a two parameter model is used to 
estimate the drag coefficient [14]. 

(2) Use of all possible analytical function matrix updates 
(elements updated analytically are those derived from the 
clock parameter differential equations). 

(3) The modified Euler integrator with two function evaluations 
per step is used. 

(4) An integration step size of six seconds is adopted. 

The results can be combined into an equivalent addition count by 
weighting the multiplications, divisions and square roots by their rela- 
tive execution times. The simulations have been performed on a CDC6600 
computer system which has the following operation times weights: 

Add 1 

Mul t 2.5 

Div 7.25 

Sq.Root 62.5 

The use of these weights with the counts of Table 4 gives the operation 
counts in Table 5 expressed in equivalent numbers of additions. Total 
value of equivalent additions in Table 5 indicates the relative effici- 
ency of the algorithms in performing covariance time updates under the 
assumptions described above. It can be observed from Table 5 that the 


square root algorithm based on Che P=UDU"^ tra-'.sform.icion in combination 
with the i propagation compares quite favorably with the conventional 
extended Kalman-Bucy filter. Note Chat if the symmetric properties of 
Eq, (15) are used, the P operation count is substantially less than either 
of Che ocher algorithms. Detailed investigations have indicated that 
there are numerical problems associated with integrating the nx(n+l)/2 
differencial equations obtained by invoking the symmetry requirements on 
P. This fact illustrates the point chat Che numerical stability is as 
important as numerical efficiency in any autonomous satellite application. 

Numerical Comparison. In addition to the operation counts described in 
the previous section, specific numerical simulations were performed on 
the algorithms to determine their actual relative performance in per- 
forming navigation computations. Tables 6 through 9 contain the results 
obtained in these simulations. The numerical results were obtained using 
Che following conditions for numerical integration of the appropriate 
differencial equations. 


(1) Variable step (2)4 Runge-KuCta with the absolute single-step 
error tolerance of 10~ , and a relative single-step tolerance 
of 10"^ 

(2) Fixed step fourth order Runge-Ku'cta with a six-second step 
size , 

(3) Fixed step, second-order Euler integrator with a three-second 
step size, and 

(4) Fixed step, secon.d-ordor Euler : r.tegrator with a stop size of 
six seconds. 


The four different integrators were used t'.i show how the relative per- 
formance changes as a function of the order and method of integration, 
the integration step size, and the method by which the integration step 
is selected. 

The data are presented in each of the tables in the following manner. The 
columns are, from leit to right; 

1. The algorithm used for covariance matrix propagation, 

2. The total measurement update time, 

3. Ihe total propagation time, 

4. The propagation time normalized by th.e lowest propagation 
time , 

5. The total computation time for time and measurement updates, 

6. The RSS of the position magnitude error, 

7. The RSS of the velocity magnitude error. 

The times are given in seconds; the position errors are expressed in meters, 
and velocity errors are given in meters per second. The algorithms are 
listed in order of increasing total enmputat ion times. 
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The resulcs indicate that the algorithm performance has a strong de- 
pendence on the characteristics of the numerical Integration algorithm. 
In general, the (4) algorithms have lower propagation times than their 
directly integrated counterparts. Tlie significant deviation from this 
rule is shown in the results of Table 9, where the direct integration 
methods shows computation time advantages. 


The other variations in performance of the algorithms are quite uneKpeeted. 
Although the computation times of all algorithm® decrease as the order of 
integration is lowered and the step size la increased, the estlaation ac- 
curacies of the algorithms increase with these same changes in integration 
characteristics. In one case, the (TO) algorithm falls to complete a 
simulation with the fourth order Integrator and a 'slx-second Integration 
step but runs successfully and competitively with the second-order inte- 
grator at the same step size. 


To rule out the possibility of computer roundoff error as a cause of this 
anomaly, tests were conducted on the Euler integrator. With step sizes 
as low as 0.5 sec., the Euler algorithm approached the results of a well 
validated high order raultistep integration code [16]. These results ruled 
out the liklihood of roundoff error or coding error as the cause for the 
anomalous characteristics in the estimation accuracy. The effect Is 
believed to be caused by a complex interaction between the stability 
characteristics of the estimator and the integrator. 

The relative performance of the algorithms in Table 9 leads Co a second 
unexpected trend in the numerical results. The predicted relative per- 
formance based on the operation counts in Table 5 does no', agree with the 
actual results obtained in Table 9. Specifically, the i..'..rerical perfor- 
mance of the 06 algorithm is better than the relative pe--’ormance pre- 
dicted by the operation counts while the W is not as good. The cause for 
this discrepancy is thought to be coding overhead not accounted for in the 
operation count or the possible parallel multiplication capability of the 
CDC6600. This question requires further investigation. 

6. Conclusions 

Based on the results given in the previous section, several general con- 
clusions can be drawn. First, the estimation algorithm performance has 
a strong dependence on the order and method used for integrating the 
differential equations involved in propagating the state estimate and 
the state estimate covariance matrix between observation epochs. This 
dependence has not been fully understood and requires further considera- 
tion, Based on the results presented here,, the square root method based 
on the UD transformation in combination wj'^-h the state transition matrix 
propagation approach appears to be the best overall square root method. 
However, for the Euler integrator using a six-second integration step, 
the best overall results were obtained with the UD algorithm. Finally, 
the single-step Euler numerical integration algorithm yields a more ac- 
curate and stable estimate than the fourth order Runge-Kutta algorithms. 

In relating the results presented in this paper to the microprocessor 
environment, several factors should be remembered. Of greatest importance 
is the fact that the timing comparisons have been obtained under conditions 




which are different than those existing for on-board computer imple- 
mentations. Table 3 and 4 are of most significance for the on-board 
application in which a high-level language, such as FORTRAN, will proba- 
bly not be used. The overhead associated with FORTRAN is the most 
probable cause of the discrepancies between the performance comparisons 
in terms of operation count versus execution times. The on-board imple- 
mentation will be in assembly language and should approach the operation 
count performance given in Tables 3 and 4, although other factors must 
be considered also. 
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Table 1. Operation Counts for ($) Algorithms 
Broken Down by Frequency. 


Algorithm 

Adds 

Mults 

Dlvs_ j 

EKF($) per update 

2n^+ 2m 

2n^+m+l 

1 

per step* 

3n^ 

2n’ 

- 

per fcn.eval. 

n* 

n' 

- 

UDU(^>) per update 

1 ,5n^+.5n^+n 
+5m+n^m 

1.5n^+2n^+.5n 
+m+ ( n ‘ -ha ) m+-l 

n-1 

per step 

3n^ 

2n’ 

- 

per fcn.eval. j 

n^ 

n^ 

- 

• i 

CARL($) per update ! 

1 . 5n^+n^+n^m 
+m 

1.5n^-ha^-.5n 

-hn^m-hn 

n-l-hn 

per step i 

3n^ 

■*» 

2n‘ 

- 

per fcn.eval. 

n^ 

i 

1 

- 

POTT(<*>) per update 

2n^+n^in+m 

2n^-f . 5n"-. 5n 
2 i 

-ha m-hn 

n-l-hn 

per step 

3n’ 

2n== 

- 

per fcn.eval. 

n^ 

n' 



Sq.Rts, 


*per step calculations assume second-order Euler integrator with 
two function evaluations. 


Table 2. Operation Counts for (P) Algorithms 
Broken Down by Frequency. 



Algorithm 

Adds 

Mults 

• 

P per step 

3 (n^-Fn) /2 

(n^+n) 

per fen. eval. 

n ‘-ha^-hn^ 


• • 

UD per step 

3 (n^-hn) /2 

n‘+n 

per fen. eval. 

n -F4n -F2n-3 

n -l-3n -ha-3 

• 

W per step 

3 (n^-ha) /2 

n 

n'+n 

per fen. eval. 

n^-fn" 

•7 ^ 

n" - . on'-f .on 


Divs Sq.Rts. 



















3 , 5n*+4ti“+ti*»+iB I 3 . 5n*+3ti* 


4n®+3u^+H*i»+« 


2n®+5Q^+3n+®* 

2n *+9 . 5u* +5 • 5tt*6 
2p'+3.5n^»1.5n 


+n a+m 


*ha*is#a 


2n^+7n *+3x1-6 


2(n +n) 


Table 4. Numerical Operation Counts 
GPS Problem. 


Algorithm 

Additions 

Multiplications 

Divisions 

Sq. Roots 

j EKF(I) 

3434 

3455 

12 

- 

UDU(I) 

3880 

4035 

22 

- 

CARLSON ($) 

3837 

3824 

30 

19 

POTTE’^d) 

4442 

4429 

30 

19 

• 

P 

3230 

2792 

- 

- 

• • 

UD 

3866 

3536 

264 

- 

W 

3102 

2684 

66 

- 





















Table 5. Numerical Operation Counts (CPS Problem). Equivalent 
Additions (CDC6600) for Modified Euler, Step Sii:e*6. 


Algorithm 

Additions 

Wghtd.Mults . 

Wghtd.Dlvs. 

Wgh'cd . 
Sq.Rts. 

Z 


3434. 

8637.5 

87. 


12156.5 

UDU($) 

3880. 

10085.0 

159.5 

- 

14124.5 

CARLSON(i) 

3837. 

9560.0 

217.5 

1187.5 

14802.0 

POTTER('l>) 

4442. 

11072.5 

217.5 

1187.5 

16919.5 

• 

P 

3230. 

6980.0 

- 

- 

10210.0 

e • 

UD 

3866. 

8840. 

264. 

- 

14620.0 

. 

W 

3102. 

6710. 

478.5 

” 1 

10290.5 


Table 6. Numerical Algorithm Comparison. Variable Step (2)4 Rung 
Kucta. Error Tolerance: REL»1D“^, aBS =■ 10“® 


Aleorithm 


Meas . 

1 1 

Norm 

Total 

Update 

Update 

Time 

Update 

8.172 

67.377 

1.000 

75.549 

8.255 

83.325 

1.237 

91.580 

8.110 

86.241 

1.280 

94.351 

18.459 

87.820 

1.303 

106.279 

9.981 

97.158 

1.442 

107.139 

20.198 

98.321 

1.459 

118.519 

7.729 

119.456 

1.773 

127.185 


Accuracy 


163.4 

163.3 
163.2 

163.2 

163.4 
163.4 

163.3 





















UDU(irD) 

CiLRL(W) 


J.SOf 
f .144 
8.015 
19,361 
9.379 
19.045 


77.014 

84.362 

84.784 

95.122 



1.000 

1.179 

1.291 

1.29S 

1.456 

1,484 


73.241 

is , 161 

92.377 
104. 14S 
104. SOI 
115.978 


MS m. 

ISl WI7 

1S3.4 

,4SI 

163.3 

.452 

163.2 

.452 

163.2 

.452 

163.2 

.431 

163,4 

,452 

- ■ 

- 


Alsorithia 


EKF 




Table 8, Numerical AlgorlCha C©ia|)sra.son , 

Modified Euler, Step SIes » 3 ««c 


Total 

Update 


72.079 

,83.565 

90.739 

102.735 

103.212 

109.803 

113.930 


Accurae 


Pos . 


7.981 

8.341 

7.867 

19.110 

9.893 

6.972 

18.818 


64,688 
75.224 
82.872 
83.625 
93.319 
102.831 
95.. 112 


1.000 

1.163 

1.281 

1.293 

1.443 

1.590 

1.470 


148.5 
146.8 
146.8 
146.8 
146.8 
146.0 

146.5 



9. Numei.cal Algorithm Comparison. 
Modified Euler, Step Size 6 sec 


Algorithm 

Meas . 

Time 

Norm 

Total 

Update 

Accuracy 

Update 

Update 

Time 

RSS Poa. 

RSS Vel. 


7,246 

38.609 

1.000 

45.555 

113.1 

.431 

EKF(^) 

8.057 

42.792 

1.108 

50.849 

120.3 

.432 

Upu(UD) 

8.045 

51.131 

1.324 

59.176 

111.4 

.433 

CAllL(W) 

19.402 

42.886 

1.111 

62.288 

126.3 

.438 

UPU($) 

8.216 

61.296 

1.588 

69.508 

117.6 

.432 

CARL(I) 

19.122 

62.404 

1.616 

81.526 

117.6 

.432 

POTT(<t) 

10.186 

71.897 

1.862 

82.083 

117.6 

.432 
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Figure 1. GPS Satellite Visibility vs. Time 



Figure 2. RSS Position Error vs. Time 
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ar* wJt* ^ f«etllea£o im*ly*i« of tli* *o|«££oii 
t* Wi 


T ( *SS 


$ 8S <i> 


M&i tii -£C'«^3N»»s^i» 

• c • 5J ^ C(£) {&5 

*. ■ ■ ■ ■:■.■■■, ... 

itse* 9 «sd 0 ia {&5 «r* io^r erlessulart awi 
»lwe« fro« {S5, C{£) i« sfc*» *y«»Beric, £lh» fol« 
lotting obaecvaEloita can be mads vegartUng (85! 

There are n(n-l)/2 unknown *le«enca in C. The pro- 
ducts u5 and Od are lower triangular eraatiftg 
a(n+15/2 unknowns. Tlierefore, eh* a » s *yst*m of 
equaelona (8) have {n(n-l5/2 ■!• a(«Tl5/2] •• a « a 
uakiMuna which can be detemlned unlqual^. 

An expansion of (8) Into matrix elsnsats 
Indicate* the method of solutioa, (tn this sxpan" 
Sion Q is_aHaumed to be a dlsgoaal matrix witk 
elements q^, • djj/2, l»!.,..,n,5 
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Sotlng that ch* first tana of (35 is che transpose 
of the aeeemd tsm. th* fsllsaliig d»|lal6l««t 1* 
mad*! 


ect) 5 


ills iafisleioq, (35*’e«;tii b* owi»r«s**d «s 


K«rh reus of th# un^ fTlsr 
mttin im m tk 

»Mbsr«d ef ^ ih* H meir- 
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■SI. t. ^ 



It 


th« B” «aerix. After a'l uppr triangular rw of 

C i« etwBpuCed, th« condition tint 

(1-1, ...,0 i l« Invoked to evalualo 

tli» corrsdpottdlng lover trlansular coluan of C. 
ThoB s eolunrn of the lover triangular •Issnertt'^ of 
M ran fe« ovoluated. Once the M snairlx elrsanitH 
are datermlnad, the next row of the upper triangu- 
lar C eloffients la cottpuled «n 1h a column of 
0 and B elemeata. Thia proceaa ts repeated until 

all ti and 5 values are determined. The U and B 
elements are determined in the following manner. 
From (63 and (7) define 


£ K + T - UB + 


C<5 


.i, " 


<-’) f i, ' A H,/,k - i T f- 
k-1 ‘ J k-H-l 


1-1 , . . . ,n ; J-H-l, . . . ,n 


1 - ! 


^l ’ ‘’ll - Vik 
1 


' -Si 


( 10 ) 


1-1 , . . . ,n : J-1 I-l 


The enponaion of K In summation notation r.lves 
n , n U,, 3, ^ 

- 1 SAt + I 
Ij ii '‘j UO ^ 


tM ci., . 1-1,. ...n 


(6) Bjj (M. J 


i-l, , . . ,n ; J*1 , , . . ,1 


( 11 ) 


1 -1 , . . . ,n ; j - 1 , , , . , 1-1 


But, since D It. diagonal, (11) becomes 

M *r T — U ci + 

"ij 'ij ij jj 2 

i-l, . . . ,n ; 1-1 , .... 1 


( 12) 


Kor 1-J , HO and D ' 1. Therefore, (12) 

becomes 


Si ■ ‘At ‘ ' " 


(M) 


lor 1 J , (12) is real r. ini', oil to I'lUalii the 

differential equation 

t 0. d 

0 » (M + T IJ-J-lwj 

IJ ^ Ij Ij 2 

i-l,...,n ; J-1 1-1 tl4) 


Equations (13) and (14) are the forma of the 
differencial equations to be employed In the 
derivative routine of a numerical integrator. lt\e 
elements and M^j are computed as outlined 

previously, and formalized In the following 
algorl tlim. 


TrlapRular Square-Root TropagaC Ion AlRorltlim 

Given the elements of tlie square-root state 
error covariance in lower triangular UD form, 

Q H Q/2, and A(t), the differential equations 

D^j and d^j^ can be computed as follows: 


i Imi ioiJ 2)2 N’ unerie.il Resul ts 

n,e iwo ri-rtinds for time pr^;;' iprit lull uf the 
squ.ir u - ri ut . i.v.i r i .trier matrix In tli.’ GDC .ilgorilhin 
were evalnated to determine the relative crjmputa- 
tli'nal spued ard in'egration accuracy. The ce.st 
prehliT. rfusi r; w is i planar Kepli'rl.in orbit .it an 
alt i Cud. ■ . t SU ',1 tn ab.ivc the oart-h. This problem 
w.is ulu s. :, ;.>r •hi. f .llowir.q reasons: 

1. lliu , 1 ,' . ■-•pltrily .illows fur qiiuk. inple- 

nri'nt.sti nr and 'rpiiT ,i t i r'n 

2 . 1 'll 1 ■ I l it, I ■ ,n ana 1 V t i I .1 1 so 1 u t ! . ri p I ves a 

sl.ind.ud ter -I’lsiiriin’. tbe aiu'urauy of the 
nitriui 1 . 1 1 !nt I .'.t ,it I . 

1. Thu or! it .■!). ..un .ipproy. Im.ites those of [iroposed 
operatloml s.itullltes (such as lAND.tAT-H) which 
will use -qo.ire-rnot algorithms In onboard navigation 
coraput urs . 

lliu tw. r. thuds wurr ev.alu.itcd by Integrating the 
st.itu ami square-root uovarianco rr.aCrl.s equ.it lotus 
for onu r.'vo I u t i on , * Tos i t i en-ve 1 oc 1 ty ur.nponents , 
statu I O..M r i am e , and lntepr.it Ion times were 
tabul.ited at Che end of the revolution as well as 
.It 1/4, 1/2, .inri 3/4 points In tbe orbit. The 
Inlegrat ion was purformod with a Kunge-K.ulCa algo- 
rithm of fourth order (RK2(4)). The time compari- 
sons vure rznde with v.irl.iMe steps Izc integration 
at relative urr ,r tnl.sranecs of lO”’, 10~' , .tnd 
If' ‘‘ . ('."iparl ons were also m.idi> for three 
fl.xed stop size.i. 

The standard ''or measuring the accuracy of 
the mimerlcal Integration of the covarl.ince matrix 
was generated by Integrating the state and square- 
root I'ov.it ian. |. uquallona at a tight tolerance 
with a high-order multi-step integrator. Covariance 
.irrur.icy w.is rroa.s.ired .is tho relative error in the 
rtignitudo of the diagonal position and vel.jcity 
rovari.incu elements. The error w.as measured 
rel.itive to the st.ind.ird solution described at the 
be ginn ing of the par.igr.iph. 

* Kumerichi comp ir i sons were porfortned un .i CDC 6600/ 
6400 uomputer ii'iing single precision -i r i l iimeC i c . 

I nt egr.i t i.'ii I , "I s .'.nil .iccur.ii: i es will v.ii/ with 
dif f ur. t t m.u li i nus . 
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1 thriiigh s, the (▼) time upd.ile algeriirr. is 

superior to the (0,D) n'.ethod i*'. enr.putaii.n 
eieney ii>r tills particular tes: pio!?len. Ipei ii'i- 
e.ally, the tollc'winR trends in relative looilthm 
rfor::.an. e hvivc* been dedin fd t I or ; the ^*■^u]ts: 

1. Fer all measwrerit?ru rates tested, the 

(J) al/'.orlthm operates faster for a ylven inte>:ra- 
l ion tolerance. Tlie differenm- in evanpul at [on 
speed incrense.s nonlinearly as i nr .»p,ra t i on 
tolerance becomes tighter. 

2. For the fixed step si.'.es an i i<o thi- 


nusasur . inent rat es t es t i‘< 
fas ter Ini ep.ra t 1 <’n time 


t hi- ( 1 I a i r ),•• 

)d h i >,h*-r O ' ai 1 •. 
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an ' • * ! i * ' e I V , h o • 
M r-.-Miod. th*' mort 


t r rei o t h*- r f lo t ■ >! ; 
should he . n*i t de T •' d 
r-n<-y St u*l i •' . 

o of the il , 1) ,» opiate 


e<Hiations in.i e/i ,it .‘d in tb.e (») algorithm. Fio 
the test pinhie-, vh»‘re the numbers of .'^na- 

tion*. to he' ir.t .*>’,rat ed. ,jre 10 and Ih, ro-pe. " i v**1 v . 
1‘erivMfive 'utin*' -.jlls still r**.piire less n.mpu- 

tativ.’.i t ir i’ 1 1 ih,* { i netloid. Perhaps a lore,. : 

state vector w'o.ld he necessary for the (L',0) to 
show a time s.nvinps due to its smaller integration 
vector size. This result will be hiRhly prohlen- 
dopendent, hut tests with larger state vectors 
woulti he en I i gh t I'n 1 nr . 

d. Th.*- (1 . 1 1 g': i t hm presented hero appr«>xl- 

nates (h.- process nois. with an analytic trapezoid 
rul.' I n f r.i • i . . Thi-. appr ox Ima t i on ren*'*v<«s an 

. id : , ■ ' 1 ! r * ■ tr.i;-; the (I> .ilgorithn. ''h.- (U.P) 


aalshiod «eee««»C» f«P pre»fe?»ii nolas **«ctly. In th» 
•inilatlotui, tb« preescs aaiea apprexlsMtlon had 
H0 dlK^ftrcAbls iCf*ct to e»v«ri«Jic» aceuracy. Pra- 
diftClTO tntanftla v«ra siMlli IncslcSv-aly, eh« 
preewa mi«« «pp<rTOiMtii»a »1U stos* 

fats M ^ 'HMMtttiMRns rag* d«er*i«Mi »t if 
elMtv* i$ a l«m pcadin&iTO latar^l »ieh m mwiwr*^ 
wmMt th* di^M of Mseume^ iggradatloa that wrold 
la an aafsaX fllcsrit^ ptobXan with Imt 
ptadia^toa la^aryal* km boon uMrasaad twra, 
bae Bfflflti stwdy. Xf th« «pproxl«#tfon w«r« not 
ta glv« tb« diwtrad aacuraey, an n x a qaadratBP* 

vwuld havt te bt addad to th« (I) al|oi?ltJ». Tha 


aiisoi i u<d lni-ri>n!ie In computatlomil hurdt-n might 

Improve the rompat Itlvenaaa of tha (U.D) updaca, 

1. Por orblf prsblaa* with highar order 
forcing funrriTOa ttoi Chat taacadt tha aecapcabla 
stapalaa wiU bacaM tmUar aa tha forcing fuoetlTO 
btcentaa lata Maeetbt Tbi' pareaaeaga reaction la 
atap a fit to accTOsiata tka anra irragular forelai 
fum-tlon nay net ba tlta aaaa for tha t’tra algorlthoiB, 
changing thalr ralatlva af flcianclaa . This prohlat- 
la rouplpd with that of tha affect of a>aaaura«acnt 
rata on atap alaa. ?ha ceablaatton of tha two 
factora datemisaa allawablt atap alae and, thara- 
fore, Intagratlcm apaad. 
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Integration Tiwaa 
25-8econc 

& Covariance 
Data Rate 

Error 



Integration Time 4 
5-Second 

Covariance 
Data Race 

Error 




T 0 L E 

R A H C E 





TOLER 

A H C E 


Rav 


10" ■« 

10‘‘ 

10‘* 

Plxed 

Rev 


30''« 

10"‘ 

10”* 

Eixad 

1 1/4 

a a 

0,D 

58858 

6240 

955 

- 

1/4 

i'.i) 

61281 

7923 

3812 

3688 


i 

9548 

1246 

751 

- 



11237 

3619 

3599 

3410 

1 

1 1/2 

9 * 

U,D 

73211 

7698 

1726 

- 

1/7 


76846 

11728 

7 576 

7409 

1 

i 

19063 

2468 

1496 

- 


: 

22612 

7251 

7203 

6860 

3/4 

U.D 

106970 

11535 

2495 

- 

3/4 

ii.D 

112596 

16886 

11404 

11103 



28590 

3671 

2229 

- 



34125 

10871 

10836 

10327 

1 

• a 

U,D 

147757 

15503 

3260 

- 

1 

U.h 

153250 

22182 

15224 

14810 



38078 

4895 

2942 

- 



45457 

14487 

14428 

13668 


COVARIANCE ERROR (1 REVOUITION) 


COVARIANCE ERROR (1 REVOIXTION) 
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M n 0 

9.74xl0“‘ 

5.14x10“’ 
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V h f' 

7.25'10'' 

2 . 73 - 10 "’ 

4. 29-10"’ 
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2.32>^10“‘ ‘ 

9.27xlo“* 
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Vel 

' 2.32'10‘‘' 

7.64x10“* 

2.73x10"’ 

4.29x10"' 

Pos 

1 1. 62x10“* 

1.61x10"* 

1.51xl0“* 

_ 

Pcs 

1 

0 

Cd 

xC 

6.25xl0“* 

6.25x10"' 

6.25x10"' 

Vel 

^ 1.46X10“* 

1.46xlo“* 

1.36x10“* 

- 
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' 5.79-10"’ 

5.79'10“* 
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Integration Times & Covariance Error 
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Integration Times & Covariance Error 
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15439 
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Baiad on tha numerical ri'bulta obtained for 
the example problem cunridered In tlie Inveat igatf on, 
Ic la concluded chat: 

— the (i) algorithm Is the more efficient In 
terms of the Integration time required to 
achieve a specified computju Ion accuracy 

• • 

— the (U,D) method can achieve the hlgheat 
computation accuracy but with a lieavy penalty 
In Integration time 




— the (U,D) algorithm has an advantage In terma 
of core storage, as It requires an integration 
vector of n(n+l)/2 elements, as opposed to the 

e 

n '' n elements required ,u the (v) method 


I'pdate sqii.- rr-roi't . ovarlancc 


- 


— the disparity In etflclemy hi'lwtvn the two 
niclhods becomes less as tin- sireul.itid measiirc- 
r.u'iit rale is Increased. 

Tlie perlurm-ince of Che 1U,D) prop.igation algorithm 
lor this test problem was not what ti.id been desired. 
Yet, results Indicate there are at le.ist two situ- 
ations wliere the (U,D) metlunl naiy efler ImgroM ii 
el f 1 clencl es. Those are In sv.stenis wltii rneasore- 
ment i.ites high enough to ci'iisuiln Int eg,: a t 1 on 
stepsii-e, and during larg.e predh tion intervals 
without measurements, where .in.ilviii- a;>pr ox ina t i an 
of tlie process noise may not be .'ulequate. These are 
the areas where future research will continue. 
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AJ'l'KNDIX 

The measurement update algorithm for the IT'!' 
I aet oil ration [7] has the following, form. '’sing 
the observation Yj^^^ j • ^ ealeul.ite; 


‘‘kdl " ^'■’^^^k+l’‘^k+l^'^^'\+l' 
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I'l - 11 + I H U 
^ ^ k-1+1 




1 • 1 ♦ n 


n+I k+1 

iiel.se) and calculate 


(where Is the measut eraai; 
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